#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stdint.h>
#include <time.h>
#include <string.h>

// 大整数类型（无符号64位，支持大范围因式分解）
typedef uint64_t bigint;

/**
 * 最大公约数（欧几里得算法）：对应Python的gcd函数
 */
bigint gcd(bigint m, bigint n) {
    if (m < n) {
        bigint temp = m;
        m = n;
        n = temp;
    }
    while (m % n != 0) {
        bigint temp = m;
        m = n;
        n = temp % n;
    }
    return n;
}

/**
 * 最小质因子筛（MPF Sieve）：对应Python的sieve函数
 * @param n: 筛法上限
 * @param prime_count: 输出参数，素数个数
 * @return: 素数数组（动态分配，需调用者释放），失败返回NULL
 */
bigint* sieve(bigint n, int* prime_count) {
    if (n < 2) {
        *prime_count = 0;
        return NULL;
    }

    // 初始化最小质因子数组
    bigint* mpfs = (bigint*)calloc(n + 1, sizeof(bigint));
    if (!mpfs) {
        fprintf(stderr, "sieve: 内存分配失败！\n");
        *prime_count = 0;
        return NULL;
    }

    // 动态数组存储素数（初始容量4，翻倍扩容）
    int pn_cap = 4;
    bigint* pn = (bigint*)malloc(pn_cap * sizeof(bigint));
    if (!pn) {
        fprintf(stderr, "sieve: 内存分配失败！\n");
        free(mpfs);
        *prime_count = 0;
        return NULL;
    }
    int pn_len = 0;

    // 核心筛法逻辑
    for (bigint i = 2; i <= n; i++) {
        bigint mpf = mpfs[i];
        if (mpf == 0) {
            // 扩容检查
            if (pn_len >= pn_cap) {
                pn_cap *= 2;
                bigint* temp = (bigint*)realloc(pn, pn_cap * sizeof(bigint));
                if (!temp) {
                    fprintf(stderr, "sieve: 内存分配失败！\n");
                    free(mpfs);
                    free(pn);
                    *prime_count = 0;
                    return NULL;
                }
                pn = temp;
            }
            pn[pn_len++] = i;
            mpf = i;
        }

        // 遍历已找到的素数
        for (int p_idx = 0; p_idx < pn_len; p_idx++) {
            bigint p = pn[p_idx];
            if (p > mpf) break;
            bigint x = p * i;
            if (x > n) break;
            mpfs[x] = p;
        }
    }

    // 释放最小质因子数组，调整素数数组大小
    free(mpfs);
    bigint* primes = (bigint*)realloc(pn, pn_len * sizeof(bigint));
    if (!primes) primes = pn;

    *prime_count = pn_len;
    return primes;
}

/**
 * Pollard Rho算法：寻找n的一个非平凡因子
 * @param n: 待分解的奇数（n>2）
 * @param p: 输出参数，找到的因子
 * @param q: 输出参数，n/p
 * @return: 成功返回0，失败返回-1
 */
int pollard(bigint n, bigint* p, bigint* q) {
    bigint b = 100000; // 初始b=10^5
    srand((unsigned int)time(NULL)); // 初始化随机数种子

    while (b <= 1000000) { // 最大b=10^6
        // 生成随机a：2 <= a <= n-1
        bigint a = 2 + (rand() % (n - 2));
        bigint g = gcd(a, n);
        if (g > 1) {
            *p = g;
            *q = n / g;
            return 0;
        }

        // 生成sieve(b)的素数列表
        int prime_count = 0;
        bigint* primes = sieve(b, &prime_count);
        if (!primes && prime_count != 0) {
            return -1;
        }

        // 遍历每个素数p
        for (int i = 0; i < prime_count; i++) {
            bigint p_prime = primes[i];
            if (p_prime >= b) continue;

            // 计算p_power = p^k（最大的k满足p^k <= b）
            bigint p_power = 1;
            while (p_power * p_prime <= b) {
                p_power *= p_prime;
            }

            // 计算gcd(p_power - 1, n)
            g = gcd(p_power - 1, n);
            if (g > 1 && g < n) {
                *p = g;
                *q = n / g;
                free(primes);
                return 0;
            }
        }

        free(primes);
        b <<= 1; // b *= 2
    }

    // 未找到非平凡因子，返回1和n
    *p = 1;
    *q = n;
    return 0;
}

/**
 * 因式分解主函数：对应Python的factor函数
 * @param n: 待分解的正整数（n>=2）
 * @param factor_count: 输出参数，因子个数
 * @return: 因子数组（升序排列，动态分配，需调用者释放），失败返回NULL
 */
bigint* factor(bigint n, int* factor_count) {
    *factor_count = 0;
    if (n < 2) {
        fprintf(stderr, "factor: n必须大于等于2！\n");
        return NULL;
    }

    // 用栈模拟Python的列表栈（最大深度32，足够分解64位整数）
    bigint stack[32];
    int stack_top = 0;
    stack[stack_top++] = n;

    // 动态数组存储因子（初始容量8，翻倍扩容）
    int factor_cap = 8;
    bigint* factors = (bigint*)malloc(factor_cap * sizeof(bigint));
    if (!factors) {
        fprintf(stderr, "factor: 内存分配失败！\n");
        return NULL;
    }

    // 栈循环分解
    while (stack_top > 0) {
        bigint x = stack[--stack_top];

        if (x == 2) {
            // 2是素数，插入到数组开头（保持升序）
            if (*factor_count >= factor_cap) {
                factor_cap *= 2;
                bigint* temp = (bigint*)realloc(factors, factor_cap * sizeof(bigint));
                if (!temp) {
                    fprintf(stderr, "factor: 内存分配失败！\n");
                    free(factors);
                    *factor_count = 0;
                    return NULL;
                }
                factors = temp;
            }
            // 数组元素后移，插入2到开头
            memmove(&factors[1], factors, *factor_count * sizeof(bigint));
            factors[0] = 2;
            (*factor_count)++;
            continue;
        }

        bigint p, q;
        if (x & 1) { // x是奇数，调用pollard分解
            if (pollard(x, &p, &q) != 0) {
                free(factors);
                *factor_count = 0;
                return NULL;
            }
        } else { // x是偶数，分解为2和x/2
            p = 2;
            q = x >> 1;
        }

        if (p == 1) {
            // 未找到非平凡因子，q是素数，插入开头
            if (*factor_count >= factor_cap) {
                factor_cap *= 2;
                bigint* temp = (bigint*)realloc(factors, factor_cap * sizeof(bigint));
                if (!temp) {
                    fprintf(stderr, "factor: 内存分配失败！\n");
                    free(factors);
                    *factor_count = 0;
                    return NULL;
                }
                factors = temp;
            }
            memmove(&factors[1], factors, *factor_count * sizeof(bigint));
            factors[0] = q;
            (*factor_count)++;
        } else {
            // 分解得到两个因子，入栈（栈是LIFO，先入q再入p保证分解顺序）
            stack[stack_top++] = p;
            stack[stack_top++] = q;
        }
    }

    // 调整因子数组为实际长度
    bigint* result = (bigint*)realloc(factors, *factor_count * sizeof(bigint));
    if (!result) result = factors;

    return result;
}

// 测试主函数
int main() {
    bigint n = 200;
    int factor_count = 0;
    bigint* factors = factor(n, &factor_count);

    if (factors) {
        printf("factor(%llu) = [", n);
        for (int i = 0; i < factor_count; i++) {
            printf("%llu", factors[i]);
            if (i != factor_count - 1) {
                printf(", ");
            }
        }
        printf("]\n");

        free(factors);
    }

    // 额外测试用例：分解561（卡迈克尔数）
    bigint n2 = 561;
    int factor_count2 = 0;
    bigint* factors2 = factor(n2, &factor_count2);
    if (factors2) {
        printf("factor(%llu) = [", n2);
        for (int i = 0; i < factor_count2; i++) {
            printf("%llu", factors2[i]);
            if (i != factor_count2 - 1) {
                printf(", ");
            }
        }
        printf("]\n");
        free(factors2);
    }

    return 0;
}